!     Copyright (C) 2010 Benjamin Piaud
!
!     LIMBES is free software; you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation; either version 3, or (at your option)
!     any later version.

!     LIMBES is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with LIMBES; if not, see <http://www.gnu.org/licenses/>

!     subroutine LIMBES_force computes at the node (xi,yi)
!     the forcing term :
!     ax, ay -> accelaration component (real*8)
!     SIMUL_gx,SIMUL_gy -> gravity component (real*8)
!     needed argument
!     integer*8 xi,yi node coordinate

subroutine LIMBES_force(ni,ax,ay)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: ni,xi,yi
  real*8 ::ax,ay,y,Ha,ay1,ay2
  Ha=1.0D-25

  xi=LIMBES_cell_xi(ni)
  yi=LIMBES_cell_yi(ni)
  y=(dble(yi)-one/two)*SIMUL_dx

  ax=SIMUL_gx
  if (y.LT.1.0d-7) ay1=Ha*(-one/y**4.0D+0)
  if (y.GT.(SIMUL_Ly-1.0d-8)) ay2=Ha*(one/(SIMUL_Ly-y)**4.0D+0)
  ay=SIMUL_gy ! +ay1+ay2 
  !      print*,xi,yi,ay
end subroutine LIMBES_force

